// --------pranmar.s---------------------------------------------- // 26-3-91 cmi pipelined ranmar generator :replaces uni() // 3-4-91 tested DL // // uses Marsaglia's RANMAR universal generator. // and gives EXACTLY same results as uni(). // u(k)=N(u(k-97)-u(k-33)), c=c-cd, if(c<0)c=c+cm, result=N(u(k)-c) // where N is normalisation 0 to 1: if(x<0) N(x)=x+1 else N(x)=x. // and cd=7654321/2^24, cm= 1.0-3*2^-24, // // also stint32, ldint32 float=> 32bit integer and converse // needed for interface with C-code rstart, seed_uni. // they are called by pransta.c and pranuni.c // // USE: // append int ustore[104]; to structure Uni_save uni_data; // void pransta(void); // void pranuni(void); // void pranmar(float *r, int n, int * ustore); // int * ustore; // rstart(i,j,k,l); // pransta(); /* converts uni_data to ustore */ // ustore = &uni_data.ustore[0]; // pranmar(r,n,ustore); /* r[0],.. contains "uni" vector */ // pranuni(); /* converts ustore to uni_data */ // // (pransta and pranuni are only needed once before several calls // of pranmar or uni respectively ). What they do is: // Convert structure: float u[98],c,cd,cm; int ui,uj; // to/from int ustore[104]; (iu[0..98],ic,icd,icm, ii,ij) // where u,c,cd,cm are mapped to 32 bit integers by // float(u) = integer *2^-32 truncated to lie 0 to (1.0 - 2^-24). // // ProtoType void pranmar( float *result, int n, int *ustore) // // Any integer n allowed: it is number of random numbers returned // // -----pranmar exactly implements the result of the following C code ----- // Uni_save uni_data; // float u[98],c,cd,cm; int ui,uj; // // float uni(void) // { /* ui=97, uj=33 at start */ // float luni; /* local variable for uni */ // // luni = uni_data.u[uni_data.ui] - uni_data.u[uni_data.uj]; // if (luni < 0.0) // luni += 1.0; // uni_data.u[uni_data.ui] = luni; // if (--uni_data.ui == 0) // uni_data.ui = 97; // if (--uni_data.uj == 0) // uni_data.uj = 97; // if ((uni_data.c -= uni_data.cd) < 0.0) // uni_data.c += uni_data.cm; // if ((luni -= uni_data.c) < 0.0) // luni += 1.0; // return luni; // } // // ii=uni_data.ui, so at start do iistart,...0 or ..(iistart-n+1) // do LOOP (n-iistart)/97 times // do tail remaining mod(n-iistart,97)-iistart // // 15 cycles per random number in loop. ==> 10.5 cycles // Timing n=500: 11.2 cycles at DL 3-4-91 C-code: 106 cycles. // // resptr = r16 n = r17 uptr = r18 ii = r19 ij = r20 c = r21 cd = r22 cm = r23 u = r24 r = r25 a = r26 b = r27 // equiv to lduptr lduptr = r27 // equiv to b uiptr = r28 ujptr = r29 step = r30 // r31 is pointer fr = f16 fr1 = f17 fc = f18 fc1 = f19 ff = f20 ff1 = f21 fs1 = f24 fs2 = f25 fs3 = f26 fs4 = f27 .align 4 _pranmar_:: // fortran load point ld.l 0(n), n // pointer to value _pranmar:: // C entry // uptr: u[98],c, cd, cm, ui, uj // 0 99 100 101 102 103 words (address * 4) ld.l 408(uptr),ii ld.l 412(uptr),ij ld.l 396(uptr),c ld.l 400(uptr),cd ld.l 404(uptr),cm adds -4,r0, step // bla decrements by 4 orh ha%const,r0,r31 // or could load orh, ixfr, f0 etc fld.q l%const(r31),fr // fr,fr1 fc,fc1 adds -4,resptr,resptr // for auto increment store //--------use ii,ij down to ii = 0 ready for loop--------------------- shl 2,ii,ii // * 4 for word addresses shl 2,ij,ij start: adds -1,n,n // test at start so 0 pass OK bc Exit ld.l ij(uptr),b subu c,cd,c ld.l ii(uptr),a bnc.t donexts // if c-cd < 0 CC clear so do next (addu) addu c,cm,c // skipped if CC set donexts: subu a,b,u subu u,c,r ixfr r,fr adds ii,uptr,lduptr st.l u,0(lduptr) // no reg+reg addressing available adds -4,ij,ij fsub.dd fr,fc,ff bc.t ij96s or 384,r0,ij ij96s: fmov.ds ff,fs1 adds -4,ii,ii // sets CC if ii < 0: so set 96 fst.l fs1,4(resptr)++ bnc start .align 8 d.pfadd.ss f0,f0,f0 // can't use d.fnop to initiate DIM nop // falls through here with ii = 96 // --------------Start of 97--------all in DIM------------------------- //--start with ii = 96, ij = 32 // do 1+8*4 passes ij=32,31-27, ...,3-0. ii=ij+64 // then do 16*4 passes ii=63-58,57-54,..,3-0. ij=ii+33 // finish with ii=96, ij=32 ready to be done. LOOP: d.fnop adds -97,n, n d.fnop bc tail // n < 97 left to do d.fnop or 28,r0,ij // bla step -4, last pass loop1 ij=-4 d.fnop adds 124,uptr,ujptr // 4*31; ujptr -16 at end of 4-loop d.fnop subu c,cd,c d.fnop ld.l 260(ujptr),a // cached loads d.fnop ld.l 4(ujptr),b // ij = 32 (31+1)*4 address d.fnop bnc.t donext00 // if c-cd < 0 CC clear so do next d.fnop addu c,cm,c // skipped if CC set donext00: d.fnop subu a,b,u d.fnop subu u,c,r d.fnop ixfr r,fr d.fnop bla step,ij,loop1 // sets LCC ij-4 d.fnop st.l u,260(ujptr) // no reg+reg addressing available //-----------------------Fast bit--8*4 passes--------------------------- // //--decrement ijptr by 16 per loop1. ij=31 to 0 done (ii=ij+64) loop1: d.fnop subu c,cd,c d.pfsub.dd fr,fc,f0 // cc:fr - cc:00 ld.l 256(ujptr),a // cached loads d.pfsub.dd f0,f0,f0 ld.l 0(ujptr),b d.pfsub.dd f0,f0,f0 bnc.t donext01 // if c-cd < 0 CC clear so do next d.fnop addu c,cm,c // skipped if CC set donext01: d.pfmov.ds ff,ff // allowed if precision same (D) subu a,b,u d.pfadd.ss f0,f0,f0 subu u,c,r d.pfadd.ss f0,f0,f0 ixfr r,fr d.pfadd.ss f0,f0,fs1 st.l u,256(ujptr) d.fnop subu c,cd,c d.pfsub.dd fr,fc,f0 // cc:fr - cc:00 ld.l 252(ujptr),a // cached loads d.pfsub.dd f0,f0,f0 ld.l -4(ujptr),b d.pfsub.dd f0,f0,f0 bnc.t donext02 d.fnop addu c,cm,c donext02: d.fnop fst.l fs1, 4(resptr)++ // can't combine 2 or 4 stores since d.pfmov.ds ff,ff // word alignment of 97 in LOOP. subu a,b,u d.pfadd.ss f0,f0,f0 subu u,c,r d.pfadd.ss f0,f0,f0 ixfr r,fr d.pfadd.ss f0,f0,fs2 st.l u,252(ujptr) d.fnop subu c,cd,c d.pfsub.dd fr,fc,f0 // cc:fr - cc:00 ld.l 248(ujptr),a // cached loads d.pfsub.dd f0,f0,f0 ld.l -8(ujptr),b d.pfsub.dd f0,f0,f0 bnc.t donext03 d.fnop addu c,cm,c donext03: d.fnop fst.l fs2, 4(resptr)++ d.pfmov.ds ff,ff subu a,b,u d.pfadd.ss f0,f0,f0 subu u,c,r d.pfadd.ss f0,f0,f0 ixfr r,fr d.pfadd.ss f0,f0,fs3 st.l u,248(ujptr) d.fnop subu c,cd,c d.pfsub.dd fr,fc,f0 // cc:fr - cc:00 ld.l 244(ujptr),a // cached loads d.pfsub.dd f0,f0,f0 ld.l -12(ujptr),b d.pfsub.dd f0,f0,f0 bnc.t donext04 d.fnop addu c,cm,c donext04: d.fnop fst.l fs3, 4(resptr)++ d.pfmov.ds ff,ff subu a,b,u d.pfadd.ss f0,f0,f0 subu u,c,r d.pfadd.ss f0,f0,f0 ixfr r,fr d.pfadd.ss f0,f0,fs4 st.l u,244(ujptr) d.fnop adds -16,ujptr,ujptr d.fnop bla step,ij,loop1 // ij-4 , d.fnop fst.l fs4, 4(resptr)++ // //----------------------------------------------------------------------- //--start with ii = 63, ij = 96 // do 64 passes ii=63-58,57-54, ...,3-0. d.fnop or 60,r0,ii // bla step -4, last pass loop2 ii=-4 d.fnop bla step,ii,loop2 // sets LCC ii-4 d.fnop adds 252,uptr,uiptr // 4*63; uiptr -16 at end of 4-loop //-----------------------Fast bit--16*4 passes--------------------------- loop2: // d.fnop subu c,cd,c d.pfsub.dd fr,fc,f0 // cc:fr - cc:00 ld.l 0(uiptr),a // cached loads d.pfsub.dd f0,f0,f0 ld.l 132(uiptr),b d.pfsub.dd f0,f0,f0 bnc.t donext11 d.fnop addu c,cm,c donext11: d.pfmov.ds ff,ff subu a,b,u d.pfadd.ss f0,f0,f0 subu u,c,r d.pfadd.ss f0,f0,f0 ixfr r,fr d.pfadd.ss f0,f0,fs1 st.l u, 0(uiptr) d.fnop subu c,cd,c d.pfsub.dd fr,fc,f0 // cc:fr - cc:00 ld.l -4(uiptr),a // cached loads d.pfsub.dd f0,f0,f0 ld.l 128(uiptr),b d.pfsub.dd f0,f0,f0 bnc.t donext12 d.fnop addu c,cm,c donext12: d.fnop fst.l fs1, 4(resptr)++ d.pfmov.ds ff,ff subu a,b,u d.pfadd.ss f0,f0,f0 subu u,c,r d.pfadd.ss f0,f0,f0 ixfr r,fr d.pfadd.ss f0,f0,fs2 st.l u, -4(uiptr) d.fnop subu c,cd,c d.pfsub.dd fr,fc,f0 // cc:fr - cc:00 ld.l -8(uiptr),a // cached loads d.pfsub.dd f0,f0,f0 ld.l 124(uiptr),b d.pfsub.dd f0,f0,f0 bnc.t donext13 d.fnop addu c,cm,c donext13: d.fnop fst.l fs2, 4(resptr)++ d.pfmov.ds ff,ff subu a,b,u d.pfadd.ss f0,f0,f0 subu u,c,r d.pfadd.ss f0,f0,f0 ixfr r,fr d.pfadd.ss f0,f0,fs3 st.l u, -8(uiptr) d.fnop subu c,cd,c d.pfsub.dd fr,fc,f0 // cc:fr - cc:00 ld.l -12(uiptr),a // cached loads d.pfsub.dd f0,f0,f0 ld.l 120(uiptr),b d.pfsub.dd f0,f0,f0 bnc.t donext14 d.fnop addu c,cm,c donext14: d.fnop fst.l fs3, 4(resptr)++ d.pfmov.ds ff,ff subu a,b,u d.pfadd.ss f0,f0,f0 subu u,c,r d.pfadd.ss f0,f0,f0 ixfr r,fr d.pfadd.ss f0,f0,fs4 st.l u,-12(uiptr) d.fnop adds -16,uiptr,uiptr d.fnop bla step,ii,loop2 // ii-4 , d.fnop fst.l fs4, 4(resptr)++ // // clear pipe-------------------------------------------------------- d.fsub.dd fr,fc,ff // cc:fr - cc:00 nop d.fmov.ds ff,ff br LOOP d.fnop fst.l ff, 4(resptr)++ //-------------------End of 97--------------------------------------- tail: // enter with ii,ij set 96,..0 range fnop adds 97,n,n // n left to do fnop or 128,r0,ij // next ij = 32*4 or 384,r0,ii // ii = 96*4 ; left DIM tailoop: adds -1,n,n // test at start so 0 pass OK bc Exit ld.l ij(uptr),b subu c,cd,c ld.l ii(uptr),a bnc.t donextt // if c-cd < 0 CC clear so do next (addu) addu c,cm,c // skipped if CC set donextt: subu a,b,u subu u,c,r ixfr r,fr adds ii,uptr,lduptr st.l u,0(lduptr) // no reg+reg addressing available fsub.dd fr,fc,ff adds -4,ij,ij bc.t ij96 or 384,r0,ij // ij = 96*4 if was negative (-4) ij96: fmov.ds ff,fs1 adds -4,ii,ii // can't have ii < 0 fst.l fs1,4(resptr)++ br tailoop nop Exit: // ii,ij are ready for next but 4* too big shr 2,ii,ii shr 2,ij,ij st.l ii,408(uptr) st.l ij,412(uptr) st.l c,396(uptr) bri r1 nop //------------------------------------------------------------------------ // 28-3-91 long stint32(float f) callable in C // // tested DL // takes 32 bit IEEE float range 0 to 0.999999 to unsigned integer // range 0 to 2^32-1. Sign of float disregarded. >=1.00 goes to .9999999. // // note ftrunc.ds of 2^32*f gives two's complement integer- so not OK // because of overflow trap. // // use IEEE float s,e(8bits),f(23bits)= (-1)^s* 2^(e-127)*1.f // _stint32:: fxfr f16,r20 andnoth 65408,r20,r21 // mask out 23 least significant bits // 111111111 0000000 = 511*128 xorh 128,r21,r21 // 1.f so mantissa *2^23, add 0 00000001 0.. shl 8,r21,r21 // 1.f mantissa * 2^31 // so store 1.f and * 2^(e-127) , typically 2^-1 or 2^-2. andh 32640,r20,r20 // e = exp +127, leading 0 (mask 0111111110..) bc zero // if f is zero too shr 23,r20,r20 // r20 should be 1,...126 (127 from 1.00 only) subs 126,r20,r20 // -(e-127 +1) (since 1.f so point 1 out) bc one_00 // CC set if 126-r20 is < 0. adds -32,r20,r0 // test underflow: 31 is max shift bnc zero // exponent too negative to give integer shr r20,r21,r16 // mantissa *2^31 * 2^(exponent+1) bri r1 // result in r16 (integer) nop underflow: // not used; just round 2^-32 to zero trap r0,r1 zero: mov r0,r16 bri r1 nop one_00: adds -1,r0,r16 // set 1111...11111 for 1.00 bri r1 nop //--------------------------------------------------------------------------- // 28-3-91 float ldint32(long i) is callable in C // tested DL // takes unsigned integer 0,..2^32-1 to float 0 to 0.99999 // so any input is acceptable. // // method double precision is s,e(1bits),f(52bits) as (-1)^s* 2^(e-1023)*1.f // so set e=1023+20=1043 = 1000001 0011, so (s,e) = 65,48 // then form (1.0+int*2^-52)*2^20 - 1.0*2^20 =result. (shift point 52 to 32) // _ldint32:: orh ha%const,r0,r31 // or could load orh, ixfr, f0 etc fld.q l%const(r31),f16 ixfr r16,f16 // 32 bit unsigned integer => 0 to 1 fsub.dd f16,f18,f20 // f17:f16 - f19:f18 // 2^20+i - 2^20 (2^52 to integer*2^-32) fmov.ds f20,f16 // result bri r1 nop .data const: .byte 0,0,0,0,0,0,48,65 // 2^20 as a double f17:16 .byte 0,0,0,0,0,0,48,65 // f19:18